Crystal nucleation for a model of globular proteins 
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A continuum model of globular proteins proposed by Talanquer and Oxtoby (J. Chem. Phys. 109, 
223 (1998)) is investigated numerically, with particular emphasis on the region near the metastable 
fluid-fluid coexistence curve. Classical nucleation theory is shown to be invalid not only in the 
vicinity of the metastable critical point but also close to the liquidus line. An approximate analytic 
solution is also presented for the shape and properties of the nucleating crystal droplet. 

PACS numbers: 



I. INTRODUCTION 



In recent years there has been an enormous increase in the number of proteins that can be isolated, due to the 
rapid advances in biotechnology. However, the determination of the function of these proteins has been slowed by 
• the difficulty of determining their crystal structure by standard X-ray crystallography. A major problem is that it 
O is difficult to grow good quality protein crystals. Experiments have clearly shown that this crystallization depends 
^ sensitively on the physical factors of the initial solution of proteins. An important observation was made by George 
and Wilson P], who showed that x-ray quality globular protein crystals only result when the second virial coefficient, 
of the osmotic pressure of the protein in solution lies within a narrow range. This corresponds to a rather narrow 
5^ ■ temperature window. For large positive B2, crystallization does not occur on observable time scales, whereas for large 
' negative amorphous precipitation occurs. Rosenbaum, Zamora and Zukoski then showed 3 the crystallization 
(~| of globular proteins could be explained as arising from attractive interactions whose range is small compared with 
O ' the molecule's diameter (corresponding to the narrow window of i?2. In this case the gas- fluid coexistence curve is 
, in a metastable region below the liquidus-solidus coexistence lines, terminating in a metastable critical point. In a 

seminal study, ten Wolde and Frenkel then showed ^3] in a simulation of a modified Lennard- Jones model with short 
T-H . range attractions that crystal nucleation is significantly enhanced in the vicinity of the metastable critical point, 
^ ' due to a reduction in the nucleation free energy barrier. Thus the conditions for optimal crystallization correspond 
' to being near the metastable critical point. The reduction in the nucleation barrier results from the large density 
fluctuations that exist near the metastable critical point, which also affect the structure of the critical nucleus. 
Namely, the initial step toward formation of the critical droplet is the creation of a dense, liquid-like droplet which 
, then forms a crystal nucleus in its interior. This process of first increasing the density and then forming a crystalline 
' order is the opposite of what normally occurs in the liquid-solid phase transition at high temperatures and close to 
, the liquidus line. The authors also showed that classical nucleation theory is invalid in the vicinity of the critical point. 



The idea that protein crystal nucleation could be enhanced by critical density fluctuations was subsequently 
■ reinforced in a density functional study by Talanquer and Oxtoby. They studied a phase field model with two 
' [ order parameters, the local density p{r) and a structural order parameter m(r). The phenomenological free energy 
functional was assumed to have a van der Waals fluid branch and a van der Waals-like solid branch. Their numerical 
solution of the Euler-Lagrange equations for the nucleating critical droplet and free energy barrier confirmed that 
^ , the nature of crystal nucleation changes qualitatively in the vicinity of the metastable critical point. They also 
' found that the nucleation barrier is significantly lower in this region than elsewhere in the phase diagram and that 
. ^ \ classical theory is invalid. Subsequently Sear studied a continuum model in which the structural order parameter 
' was omitted 4\. He carried out a mean field theory near the metastable critical point and showed that the liquid 
[ "tail" of the density profile of the critical droplet is described by an Ornstein-Zernike like exponential behavior, 
whose range is given by the fluid correlation length. Since this length diverges as one approaches the critical point, 
so does the critical size of the droplet and the excess number of particles within it. In a subsequent paper Sear 
extended his arguments via a scaling theory using the critical exponents for the Ising universality class, so as to be 
valid in the vicinity of the critical point. In particular, the singular part of the free energy barrier associated with 
the tail was shown to satisfy the same scaling as the order parameter associated with the metastable critical point 
(the density difference) and vanishes as one approaches Tc. Sear did not calculate the shape of the droplet in the 
core and interface region and hence was unable to calculate the value of the barrier to nucleation. 

In this paper we extend the work of Talanquer and Oxtoby, to obtain a better understanding of homogeneous crystal 
nucleation. We present additional numerical results for the model as well as an approximate analytical description of 
the critical droplet, including its crystal core and tail (the latter is based on Sear's work Q). From this we compute 
the free-energy barrier to nucleation in the region between the fluid-solid coexistence curve and the metastable fluid- 
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fluid coexistence curve. As expected, this barrier is smallest in the vicinity of the metastable critical point. We show 
that classical nucleation theory is incorrect in the region between the liquidus and solid coexistence lines, except 
presumably very close to the liquidus line. 

The outline of the paper is the following. In section 2 we provide a brief review of the model. In section 3 we 
extend the numerical analysis of for the nucleation free energy barrier along various thermodynamic paths. We 
also calculate the profile of the critical droplet as a function of the density and structural order parameter. This 
yields the excess number of molecules, N , and the number of molecules in the crystalline phase, Nc ,of the critical 
droplet. In the next section we present a brief discussion of the range of validity of classical nucleation theory and 
show that it is incorrect in the entire metastable region of the phase diagram that we explore. We also show results 
for the surface tension along different paths in the phase diagram. In section 5 we present an approximate theory for 
the shape of the critical droplet and from this obtain N and Nc- The theory accurately describes the critical droplet 
except in the interface region, and is in qualitative agreement with the numerical results for N/Nc- Section 6 contains 
a brief conclusion. Some details of the numerical analysis are given in an Appendix. 



II. MODEL 



We use a model due to Talanquer and Oxtoby ^ to describe globular protein crystallization. Their phase field 
model is based on the following grand canonical free-energy functional: 



Q [p, m] — J 



f{p,m) - pp+^KpiVpf + ]^K„,pl{Vmf 



(1) 



The free energy depends on two order parameters: the (conserved) local density p(r, t) and a (non-conserved) local 
structural order parameter that shows whether the system is in a solid or fluid phase m(r,t). Here f{p,m) is the 
Helmholtz free-energy density and is the chemical potential. Talanquer and Oxtoby P| use the van der Waals free 
energy density for the fluid branch: 

ff{p, m) = ksTp [Inp — 1 — ln(l — pb)] — ap^ + ki,Taim^ (2) 

and a corresponding van der Waals free energy for the solid phase: 

/s(p, m) = ksTp [Inp - 1 - ln(l - pb)] - (a + araml)p^ + ki,T{aim^ + 02) (3) 

Here p plays a role of a weighted coarse grained density p = p{l — a^m{2 — m)) where the term in square brackets 
implements a difference between the fluid and solid close-packing limits. For the fluid close-packing limit, pb = 1, 
while for the solid close-packing limit pb — 1. The quantity ms{p) is the equilibrium value of the order parameter 
in the solid phase and hence is a solution of the equation {dfs/d'm)m^(p) — 0. Thus in the solid close-packing limit 
{bp = 1) vTLs = 1. The parameter am in equation ijSJl has been introduced in order to change the range of the attractive 
interactions between molecules in the solid phase from that in the liquid phase (as described by the parameter a in 
equation |2Jl). 

The chemical potential in the solid and fluid phases is the first derivative of the free energy with respect to density: 
p — {df/dp)T- In order to get coexisting densities pa and pp we have to solve the equations p{pa) = pip/3) and 
^{pa) = ^{Pi3)i where uj = f — pp. Graphically this gives the well-known "common tangent" rule for coexistence. By 
repeating this for different temperatures we can obtain the entire phase diagram. 

In order to obtain a critical droplet profile we must solve the Euler-Lagrange equations with appropriate boundary 
conditions: 

6n , m 

= and -— = 
Op dm 



-Kp\7^p+^-p^0 (4) 
om 

Using the solutions of Q for p(r) and m(r) we can obtain such properties of the inhomogeneous system as the free 
energy barrier AfJ — ri{p(r), m(r)} — il{po, 0} and the surface tension cr = 1/2(/q°° iCp I Vpl dr + Kra\'^ra\ dr). 
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FIG. 1: Phase diagram for the particles with the short-range attraction interactions. The fluid-fluid coexistence curve becomes 
metastable in this case. 



III. NUMERICAL RESULTS 



In Appendix A we describe the numerical methods that were used to solving equations Q. Here we present some 
results, extending the work of 0, for the following choice of parameters: a ^ 1, b = 1, ai = 0.25, a2 = 2, — 0.3, 
Gm = 1- The phase diagram for these values is shown in Figure 1, which shows in particular a metastable fluid-fluid 
coexistence curve. The existence of the metastable critical point affects the nucleation and growth processes in the 
vicinity of this point. 

The main quantity characterizing the nucleation process is a nucleation rate /, given by 

/ = /oe-^^/'^^ (5) 

where Jq is a prefactor given by the product of dynamical and statistical parts 0- In order to understand how this 
metastable critical point affects nucleation, we calculate the free energy barrier for different thermodynamic paths. 
(Some results for this barrier were presented in @.) The barrier dependence on temperature and density in the 
metastable region between the liquidus (solubility curve) and solidus lines is shown in Figures 2 and 3. Figure 2 shows 
that as the temperature decreases, the barrier also decreases, while Figure 3 shows that as the density increases, the 
barrier decreases. This is the normal behavior of a free energy barrier, since at any point on the solubility curve the 
free energy is infinite, decreases as one moves away from it and vanishes at the solid coexistence curve. 

However, the existence of the metastable fluid-fluid coexistence curve changes the pathways of the constant free 
energy barrier and constant supersaturation lines (Fig. 4) as compared with the case in which the coexistence curve is 
not metastable. We can also see that along the lines with constant supersaturation the free energy barrier decreases 
as one approaches the critical point (Fig. 5), which is consistent with 6]. Figure 5 shows that the free energy barrier 
has a minimum near, but not at, the critical point. The location of this minimum changes from above the critical 
point for low supersaturation to below the critical point for high supersaturation. It also can be seen that the increase 
of the free energy barrier below the critical point becomes very sharp for the cases where the constant supersaturation 
lines intersect the fluid-fluid coexistence curve 0. In order to understand the behavior of the free energy barrier 
along the constant supersaturation lines, we flrst consider their shape (Fig. 4). At high temperatures these lines are 
more vertical (constant density lines) and the free energy barrier decreases with temperature, as shown in Fig. 2. 
Near the critical point the curves become almost horizontal (isothermal lines) and the barrier starts to increase (Fig. 
3) . Thus the shape of the constant supersaturation curve in some sense determines the behavior of the free energy 
barrier. To understand this in more detail, we note that the derivative of temperature with respect to density at 
constant supersaturation is 

rdT\ (^) 



J pp 

dT i 



(6) 
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FIG. 2: Dependence of the free energy barrier versus temperature at constant density. 
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FIG. 3: Dependence of the free energy barrier versus density at constant temperature. 

where fpp is evaluated at the background fluid density. This derivative vanishes at the critical point, so that the 
constant supersaturation lines become horizontal in the vicinity of the critical point. On the other hand, at large 
temperatures the denominator vanishes, so that the constant supersaturation lines become vertical. Thus we see 
that the presence of the critical point changes the behavior of the free energy barrier from decreasing as we lower 
the temperature far from the critical point (the vertical part of the A/i — const lines) to increasing as we lower the 
temperature to its critical value (the horizontal part of the A/x = const lines). Therefore somewhere in between there 
is a minimum of the free energy barrier. Thus we can conclude that the free energy barrier is relatively unaffected by 
the existence of the critical point along paths of constant temperature or constant density , whereas it plays a crucial 
role along paths of constant supersaturation 

One quantity of interest is the excess number of molecules, defined as the number of molecules in the presence of 
the droplet relative to the number of molecules in the spatially homogeneous metastable state: 



N = f {p{r)-po)dr 
Jv 



(7) 



5 




FIG. 4: The solid lines on the phase diagrams show the paths of constant supersaturation. The dashed line shows the border 
between the background fluid and background solid. If the background density is to the left of the dashed line, it corresponds 
to the fluid branch of the free energy. If it is to the right, it corresponds to the solid branch. 
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FIG. 5: Dependence of the free energy barrier versus temperature at constant supersaturation. 



We can calculate the number of molecules in the crystaUine state using m(r), since this phase field shows whether a 
point is in the solid or fluid state. Thus: 

Nc= [ m{r){p{r) - po)dr (8) 
Jv 

Mean field theory yields a divergence in the excess number of particles at the metastable critical point, whereas the 
number of particles in the crystalline state remains finite. Figure 6 shows Nc versus N for different supersaturations. 
We can check our numerical results with the nucleation theorem lol llOl : 



(9) 
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FIG. 6: The number of crystalline particles versus excess number of particles for different temperatures. 
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FIG. 7: The excess number of particles versus supersaturation (solid line). The circles axe obtained by taking the derivatives 
of the free energy barrier with respect to the supersaturation 

Thus wc can take a derivative of the free c;nergy barrier with respect to supersaturation and compare this with our 
results for the dependence of the excess number of particles on the supersaturation (Fig. 7). We see that the free 
energy barrier decreases rapidly near the critical point as a function of the supersaturation. This happens because 
the background fluid density as a function of the supersaturation becomes flat in the vicinity of the critical point. 



IV. BREAKDOWN OF CLASSICAL NUCLEATION THEORY 



Classical nucleation theory (CNT) assumes that the critical droplet is large and that its interface is very sharp. It 
predicts that the free energy barrier is 



Mid = 



3(p, - po)2A/x2 



(10) 



7 




FIG. 8: Density profiles for different temperatures. From left to right: T=1.01, T=1.25, T=1.5. The horizontal line shows the 
background fluid density. 



where A/z is /i — ^coex- However, classical nucleation theory becomes increasingly inaccurate as one approaches the 
metastable critical point This is because the critical droplet has a diffuse interface in this region, in contradiction 

with the assumption of CNT. Our goal in this section is to compare CNT with our numerical results for the free energy 
barrier and to show that CNT has at best a very limited region of validity. It is useful to divide the contributions 
to the free energy barrier into two parts, i.e. the bulk and surface contributions, and compare our numerical results 
with the predictions of CNT for each of these separately. First, define 



Af2 



hulk 



4:7T / Auj{r)r'^dr, 



(11) 



An, = 47r y i [Kp\S7p\^ + K^\Vm\^] r^dr 



so the free energy barrier can be written as 



An = An 



bulk 



An„ 



The classical approximation for this barrier is 



AuJr: 



(12) 



(13) 



(14) 



In the above Aw(r) = [f {p{r) , m{r)) — fJ,p{r)] — [/(po,0) — ppo] is the difference between the value of the grand 
canonical free energy at the points r and infinity, Aujci — [/(Ps,Ws) — pps] — [f{po,0) — MPo] and a is the surface 
tension. This division also allows us to determine which of these contributions to the free energy barrier is more 
important in the breakdown of CNT. We first note that if CNT is correct at some point {p, T) of the phase diagram, 
then the first terms of and (|14|l should be equal; likewise, the second terms also should be equal. Therefore the 
values of the radii that result from comparison of the first terms, i.e. i?i = (3Af2b„;fe/|A(jJc/|)^^'^, as well as the second 
terms, R2 = {An, / uY^'^ should be equal to each other and to the classical radius of the droplet, Rd = 2(t /\AijJci\- 

Using this observation, we have checked the accuracy of CNT close to the liquidus line at low temperatures, as well 
as close to the critical point. The point close to the critical point is T/Tc = 1.01 and ph = 0.33, and those close to the 
liquidus line are T/Tc = 1.5, pb = 0.33 and T/Tc — 1.4, p = 0.2. For all these points CNT is incorrect. (For reasons 
described at the end of the Appendix we have been unable to explore the region very close to the liquidus line, where 
CNT presumably becomes correct.) Figure 8 shows the density profiles for these cases, while Table 1 shows the values 
of the corresponding quantities. An, An^, R2, Rci and the actual radius of the droplet Ra- We define the latter 
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100.7 


12.3 


3.83 


6.11 
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6.53 


Close to critical point 


3.49 


0.28 
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1.28 
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4.09 



TABLE I: Values for comparison of our results with the classical nucleation theory. The points close to the critical point and 
liquidus line are T/Tc = 1.01 and pb = 0.33, and T/Tc = 1.5 and pb = 0.33, respectively. 
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FIG. 9: Grand canonical free energy density versus r for two different points on the phase diagram. 

as the radius of a spherical droplet which has an infinitely sharp interface and the same solid density and number of 
particles as our droplet. As can be seen from Table 1, CNT underestimates the nucleation barrier by a factor of ten as 
compared with the actual value for the model. It also underestimates the critical droplet size by a factor of two. As 
one would expect, the ratio of the actual and classical nucleation barriers decreases as we approach the liquidus line. 
However, it is important to note that even in the region which is 'close' to the liquidus line the difference between the 
actual and classical values of the barrier is still large. It seems, therefore, that CNT is accurate only very close to the 
liquidus line. 

Next we try to determine how CNT breaks down. From Table 1 we see that the radius i?2 extracted from the surface 
contribution to the free energy barrier is closer to the actual value of the droplet radius than the radius Ri extracted 
from the bulk part. In order to understand why Ri is significantly smaller than the actual radius, we examine the 
dependence of Alu on r (Fig. 9). At the solid core Au is negative and equal to the classical AlUcI = i^soHd — ^ fluid 
(dashed line in Fig. 9). As we reach the interface region of the droplet, Aoj starts to increase, becomes positive and 
reaches its maximum at the inflection point of the density profile As r increases beyond the interface region, 

Aco decays and approaches zero as we reach the background value of the fluid density. In order to obtain the bulk 
contribution to the barrier one must multiply Auj by and then integrate over r. The dependencies of Auj{r) and 
the 'incomplete' barrier AQ{r) — /J Auj{r) *r^dr on r are shown in Figures 10 and 11, respectively. The dashed lines 
in these figures are Acuci * and An /ir^ Aojci- The dotted line in Figure 11 is the value of the bulk contribution to 
the free energy barrier (which is equal to the value of the incomplete barrier as r — > oo). The value of Ri corresponds 
to the radius at which the bulk contribution to the classical free energy is equal to the bulk contribution of the actual 
free energy barrier (given by the intersection of the dotted and dashed lines in Fig. 11). 

Therefore, one can see that since the interface of the droplet is not sufficiently sharp, the radius and therefore 
the classical radius Rd , is significantly different from the actual radius, Ra^ of the droplet. This is responsible for 
the large difference between the classical and nonclassical values of the free energy barrier. The main contribution to 
the interface is its " tail" , which is given by the correlation length and hence becomes very large close to the critical 
point. At temperatures close to the critical temperature the system should be in the very narrow region near the 
fluid coexistence curve in order for the radius of the droplet to be much larger than the interface thickness, as given 
by the correlation length. 
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FIG. 10: Grand canonical free energy density times (integrand) versus r 
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FIG. 11: 'Incomplete' bulk contribution to the free energy barrier versus r. Dashed line - dependence of the bulk contribution 
to the classical free energy barrier on the classical radius of the droplet. Dotted horizontal line is the bulk contribution of the 
free energy barrier. 

V. THEORETICAL ANALYSIS 
A. Approximations for the critical droplet. 

Equations Q with appropriate boundary conditions define the saddle point of the free energy landscape in functional 
space. The boundary conditions are a) that the density and structural order parameter at infinity are given by the 
metastable values and b) their spatial derivatives are zero at the origin, i.e. p(oo) = po, m{oo) — 0, dp/ dr{Q) = and 
dm/dr(0) = 0. There is a trivial, spatially uniform solution, satisfying these boundary conditions, given by p(r) = pg, 
m{r) = 0. This uniform solution corresponds to the supercooled liquid, which is a metastable "point" in functional 
space. The critical droplet is a spatially inhomogeneous, unstable solution satisfying the same boundary conditions, 
which we find by first solving the Euler-Lagrange equations around the droplet core (r = 0). We assume that at the 
origin the density has a value close to ps and the structural order parameter has a value close to rUg ■ [l^ (In general 
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the values of the density and structural order parameter at the droplet center are not equal to their equilibrium values, 
in contrast to the assumption of classical nucleation theory.) First, we represent the solution in the following form: 



p{r) ^ ps-i 
m{r) — rris + 



Xpjr) 
r 

Xm{r) 



(15) 
(16) 



We substitute H15|l and (|16|l into rewrite these equations in terms of the x's and then linearize them to obtain 



K, 



-K„ 



^ Am 



~^ fppXp ^" fprnXrn — 
~l~ frnp\p ~\~ fmmXm — 



(17) 
(18) 



Here fpp, and /,„ are the second derivatives of f at pg and mg. Because this is an approximation near the droplet 
core, we obviously cannot use the boundary conditions at infinity. So, in addition to the zero derivatives at r = 0, 
we add two arbitrary boundary conditions at the origin. The boundary conditions for the x's therefore become the 
following: Xp(0) = 0, Xm(0) = 0, Xp(0) = and Xm(0) = Am, where Ap and Am are the (unknown) shifts in the 
values of density and order parameter from ps and rus- The values of Ap and Am cannot be equal to zero, because 
in this case the solution is spatially uniform (see Appendix A). Solving equations (17) and (18), we obtain 



sinh{q2r) sinh{qir)\ Ap ( ql . . 1i ■ \ 
Xp = A\ + 5" — sinh(qir) sinh{q2r) 

' <72 qi J qi-qi \qi 



92 



and 



where 



smh{q2r) sinh{qir)\ Am ( ql . ql ■ . 
Xm = B + 2 —smh[qir) smh(q2r) 

92 qi J qi- qf \qi 92 



^ _ fppAp + fpmAm ^ _ f„ipAp + franiAm 



(19) 



(20) 



Kp{qi-qi] 



Km{qi~qi] 



fpp , fn 



Kp K„ 



± 



fpp fn 



Kp K„ 



4f2 

J pm 

KpKrn 



(21) 



Thus we obtain the solution at the center of the droplet in terms of two arbitrary parameters, Ap and Am. In principle 
these parameters could be obtained from the exact solution by matching this form to the boundary conditions at 
infinity. 

Next we find the solution for the tail of the droplet 3 . As r — > cxo we have metastable boundary conditions, so we 
choose for the solution small deviations from the metastable values: 

p{r) =po + ^e-'"^ 



m(r) — — e ^ 
r 



(22) 



where Sp and 6,, 



are arbitrary parameters and qp = l/C^' = {Kp/fpl^y^'^ and qm = l/^£^ = (i^m//mm)^^^- The 
tail solution for m is in fact exact. We now have solutions in both the core and tail regions, but not in the interface 
region. Figure 12 shows the comparison of the numerically obtained density profile at T/Tc = 1.2 and po = 0.33 
with our approximations for the tail and the core. The tail approximation is rather good, but the core approximation 
becomes incorrect as one approaches the interface region. For comparison purposes, we used the values of Ap, Am, 
dp and dm obtained numerically. 



B. Approximation for N^/N. 



N and Nc can each be considered as the sum of three terms: the number in the core, the tail and the interface. We 
can calculate the contributions to each of these from the core and tail, using the solutions for these two regions from 
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FIG. 12: Comparison of the numerically obtained density profile (at T /Tc = 1.2 and po = 0.33) with the tail and core 
approximations (dashed lines). 



the above section. Hence we can write 



Jo ^ 



(23) 



-r dr 



(24) 



Rta 



where Rcore is the maximum value of r for which we can still use the approximation near the center of the droplet, 
and Rtaii is the minimum value of r for which the tail approximation can be considered as correct. The interface 
region is Rcore < r < Rtaii- In the core region, << Ps so we can drop the second term in (|23|l . Integrating H23|) 
and 1211 we find 



and similarly 



j^icore) ^ _ Po)F('=°'-'=) 



(ms H )[ps + po)r dr 

r r 



(25) 
(26) 

(27) 



which gives us 



nritail) _ X X 



Opt- "mt; 2 , 

— r dr 



Rtai 



Next, we assume we can neglect contributions from the interface region, i.e. 

n!:'---) + 7vj«"0 ^ Nc 

]\[{core) _|_ ]\[{tail) ^ ]\[ 



(28) 

(29) 
(30) 

(31) 
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FIG. 13: Approximation for the ratio of the excess number of particles and the number of crystalline particles versus background 
density for T/Tc = 1.01, obtained by matching core solution directly with tail solution. 



In order to check this approximation, we compare its results with the numerical results. However, to do this we must 
discuss how we choose to determine the droplet profile, in order to estimate the ratio of Nc/N. 

We first note that the LHS of (|31|l depends on two parameters, Rcore and Rtaii (for simplicity we didn't introduce 
different radii for p and m). Rcore depends on the amplitudes Ap and Am, while Rtaii depends on the amplitudes 
Sp and 5m- Therefore the LHS of (|31|l depends on four parameters. If we knew the interface solution analytically, 
we would be able to determine these four parameters by matching these solutions. As the approximation we could 
match the core and the tail solutions (so Rcore = Rtaii)- This is the easiest to do but the least accurate. By setting 
Rcore = Rtaii = Rmatch we must take iuto account that the core solution is obtained for the solid branch and the tail 
solution is found for the fluid branch of the free energy. Therefore Rmatch is a point where the solid branch of the 
free energy intersects the fluid branch. 

Figure 13 shows the dependence of the ratio of the excess number of particles and number of crystalline particles 
on the background densities for the case in which we match the core solution directly with the tail solution. As can 
be seen, this method gives the best estimate for this ratio near the critical point, since the tail plays a more important 
role in that region and the approximation for the tail of p is more accurate there. However, the actual values of N 
and Nc are each underestimated by a factor of ten . 

As an alternative method we use the numerically determined Ap and Am to determine the core solution. Then 
we match the tail solution with this core solution at the point of intersection of the fluid and solid free energies. In 
this method the derivatives of p and m with respect to r are not equal at the fitting point (Fig. 14). As can be seen 
from Figures 15 and 16 this approximation is qualitatively good close to the critical point, but becomes worse as T/Tc 
becomes larger. We still underestimate N and Nc, but not by as large a factor as in the above case. 

However, close to the metastable critical point 7V(*''*') (and therefore N) diverges and becomes the leading term 
in the ratio (|31|l . while Nc stays finite. Thus we can obtain the asymptotic behavior of the ratio (|31|l close to the 
metastable critical point. In this region we can rewrite H31|l as 

^ « oc ^ oc f(/) (32) 

Therefore close to the metastable critical point this ratio vanishes linearly with r and quadratically with Apc'- 

'^'^ «(r+f(Ap.)^) (33) 

metastable \ / 



TV 



where t = {T — Tc)/Tc and Apc = po — Pc- Another asymptotic limit is the approach to the liquidus line. In this case 
the radius of the solid core diverges and therefore N'^'^"^'^^ and Nc"^"^*^^ become leading terms in (27) so at the liquidus 
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FIG. 14: Comparison of the numerically obtained density profile (at T/Tc = 1.2 and po = 0.33) with the solution obtained by 
matching core solution (with parameters Am and Ap obtained numerically) with tail solution, without matching the derivatives. 




FIG. 15: Approximation for the ratio of the excess number of particles and the number of crystalline particles versus background 
density for T/Tc = 1.01 obtained by matching the core solution (with parameters Am and Ap obtained numerically) with the 
tail solution, without matching the derivatives. 



line the ratio becomes 

= (34) 



N 



solubility 



We can see that there are two cases where the number of molecules in the critical droplet diverges, namely close 
to the critical point and close to the liquidus line, respectively. However, the number of molecules in the solid core 
of the droplet diverges only in the latter case. Thus close to the critical point most of the molecules are in the tail 
of the droplet, due to the divergence of the fluid correlation length, as first pointed out by Sear At the liquidus 
line most of the molecules are in the solid core, whose radius diverges. Thus instead of one length scale (the radius of 
the droplet) we have two: the radius of the solid core and the fluid correlation length. The solid core of the globular 
protein crystal droplet can be observed experimentally , , where it has been shown that close to the metastable 
critical point the core consists of only a few molecules( about 4-10) . However, the existence of the long tail in the 
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FIG. 16: Approximation for the ratio of the excess number of particles and the number of crystalhne particles versus background 
density for T/Tc = 1.05 obtained by matching the core solution (with parameters Am and Ap obtained numerically) with the 
tail solution, without matching the derivatives. 



vicinity of the metastable critical point has not been directly observed. 



In this paper we have extended the numerical results |a] to obtain a better understanding of nucleation for their 
model. In particular, we have calculated the density and structure order parameter profiles of the critical nucleating 
droplet for different temperatures. This solution of the saddle point equations then yields the free energy barrier to 
nucleation, which we have shown for a variety of different paths in the phase diagram. We also show in detail how 
classical nucleation theory is invalid everywhere in the region we studied, except presumably very close to the liquidus 
line. An approximate theory for the shape and properties of the critical droplet is developed that is in reasonable 
agreement with the numerical results, although there is clearly a need for a better description of the interface region 
that separates the core and tail of the droplet. Finally, we note that experiments show the existence of a gel state 
for globular proteins, in part of the metastable region of the phase diagram E6\. As a consequence, our discussion 
of the free energy barrier to nucleation for globular proteins would only apply to the region of the phase diagram in 
which there is no gelation, which is (roughly speaking) for densities p < Pc- The model studied here does not exhibit 
gelation. 



We wish to thank D. Oxtoby, V. Talanquer and P. Vekilov for helpful discussions. This work was supported by a 
grant from the National Science Foundation, DMR-0302598. 



Solving the equations I0J is a nonlinear two-boundary problem. One of the standard methods of solving this type of 
equation is the shooting method. Because is it very difficult to use the conditions at both boundaries simultaneously, 
we implement the boundary conditions on one side, and choose the number of free parameters equal to the number 
of boundary conditions on the other side. Then we can solve the resulting initial value problem. Depending on the 
final values (on the second boundary), we can decide how to change these free parameters. In our case we have two 
boundary conditions for r = 0, which we can use directly, and two boundary conditions at r = cxd, which we cannot 
use directly. However, we can use Ap and Am as free parameters. 



VI. CONCLUSION 
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VIII. APPENDIX A. NUMERICAL SOLUTION OF SADDLE POINT EQUATIONS. 
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We still have several difficulties. First of all, the second boundary is at infinity, so we can't explicitly get the 
values on this boundary. We could use some large finite value instead of infinity and avoid this problem, as long as 
the solution doesn't diverge at infinity for an incorrect choice of the parameters Ap and Am. Unfortunately it does 
diverge, so the linear shooting method doesn't work in this simple way. The reasons for the divergence are the error 
involved with the numerical method and the absence of solutions for most choices of Ap and Am. The error involved 
with the numerical method increases as we approach a cusp in free energy landscape; this is one of the reasons for the 
divergence. By varying the parameters Ap and Am we can push the beginning of this divergence to larger distances, 
so we can obtain part of the solution (including part of the tail) before the divergence occurs. However, this gives large 
errors in estimating such quantities as the free energy barrier, surface tension, etc. We could decrease the numerical 
error of the shooting method in the interface region by using Runge-Kutta methods with adaptive step size, but if we 
desire reasonably small values of the errors, the step size becomes too small for simulations [l7| . 

To avoid these difhculties we use a more advanced shooting technique - shooting with a fitting point. The fitting 
point was chosen at r = Rcore where the solid and fluid branches of the free energy intersect. The value of Rcore 
depends on the parameters Ap and Am. Next, we know the exact solution for the fluid part of m(r) (|22|l . where Sp is 
a parameter. By matching the fluid and solid parts of m(r) and their derivatives at r = Rcore we determine Am for 
a flxed values of Ap and 6p. The next step is to obtain Ap and 6p by matching the values and derivatives of the fluid 
and solid parts of p{r). Note that the value of Am depends on the value of Ap, so when we shoot with a particular 
value of Ap we obtain a value for Am each time. This shooting with a fitting point minimizes the numerical error as 
compared with standard shotting method. 

There still is one remaining problem, even if we use shooting with a fitting point. Namely, as we approach the fiuid 
coexistence curve, the values of Ap and Am become very small, so for shooting and fitting we need to vary them 
very slightly. Thus it becomes beyond the numerical limit of precision of the calculations to obtain a matching at the 
fitting point as we change these parameters. This so far prevented us from exploring the free energy barrier as one 
approaches the coexistence curve and hence we have not been able to determine where the classical theory becomes 
valid. 
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